%% Make Table S1

% Run after running the first section of 'calculate_flagella_basis_modes.m'

%%
try
for i = new_n_allruns%1:length(allruns)
    numFrame = size(allruns(i).tail_final,2)-7;
    xy1 = allruns(i).tail_final{4}(1,1:2);
    xy2 = allruns(i).tail_final{end-4}(1,1:2);
    
    new_vel(i) = sqrt(sum((xy2-xy1).^2))*um2px*params(1).time_res/numFrame; 
end
catch
end

%% Velocity (Swimming Speed)

vel = [];
for i = new_n_allruns%1:size(allruns,2)
    %     weight(i) = size(allruns(i).head,1);
    %     vel(i) = allruns(i).vel;
    vel = [vel;allruns(i).newVel];
    
    
end

% MeanVel = sum(weight.*vel)/(sum(weight))
MeanVel = nanmean(vel);
MedianVel = nanmedian(vel)
StdVel = nanstd(vel)
% length(vel)

%% Flagellar length

idx= 0; len = 0;
for i = 1:size(allruns,2)
    for j = 1:size(allruns(i).tail,1)
        
        
        %         tailTemp = allruns(i).tail{j}(:,1:2);
        %         dlength = diff(tailTemp);
        %         dlengthMag = dlength(:,1).^2+dlength(:,2).^2;
        
        tailFlag = (allruns(i).tail{j, 1}(:,3) > 0.001 | allruns(i).tail{j, 1}(:,3)==0); tailFlag = tailFlag(2:end);
        if sum(tailFlag) == 0
            len(idx+1) = NaN;
            
            idx = idx+1;
            
        else
            tailFlag = [tailFlag,tailFlag,tailFlag];
            
            % Hack Job. This works, don't worry about it
            CC = bwconncomp(tailFlag);
            tailFlag(:)=0;
            numPixels = cellfun(@numel,CC.PixelIdxList);
            [biggest,nn] = max(numPixels);
            tailFlag(CC.PixelIdxList{nn}) = 1;
            tailFlag = tailFlag(:,2);
            nn = find(tailFlag,1,'first');
            tailFlag(1:nn) = 1;
            
            
            inc_length = diff(allruns(i).tail{j, 1}(:,1:2));
            temp_length = sqrt(sum(inc_length.^2,2));
            len(idx+1) = sum(temp_length(tailFlag));
            temptailLength(j) = sum(temp_length(tailFlag));
            idx = idx+1;
        end
    end
    tailLength(i) = nanmean(temptailLength(:))*um2px;
    temptailLength = [];
end



avg_length = nansum(len)*um2px/idx
std_length = nanstd(len*um2px)

%% Cumulative 2nd mode energy

%Table
for i=size(S1,1):-1:1
    if sum(S1(i,:)) == 0
        S1(i,:) = [];
    end
end
cumEnergy = nanmean(S1,1);
cumStdev = nanstd(S1,[],1);
cumEnergy(2)
cumStdev(2)


%% Wavelength
    y=[];
for nidx=new_n_allruns
    y_ = zeros(30*size(cell2mat(allruns(nidx).curvature),2),1); 
    for i = 1:size(cell2mat(allruns(nidx).curvature),2);
        y_((30*(i-1))+(1:30)) = allruns(nidx).curvature{i};
        %y = [y;allruns(nidx).curvature{i}];        
    end
    y = [y;y_];
%     y = cell2mat(allruns(nidx).curvature);
    %Fs =  params(1).time_res; %Hz
    Fs = 30;
    
   
end
 NFFT = length(y);
    Y = fft(y,NFFT);
    Y = fftshift(Y);
    dF = Fs/NFFT;                      % hertz
    F = -Fs/2:dF:Fs/2-dF;           % hertz
    magnitudeY = abs(Y);
    phaseY = unwrap(angle(Y));
    %helperFrequencyAnalysisPlot1(F,magnitudeY,phaseY,NFFT)
    plot(F, magnitudeY); xlim([0 max(F)])
    
    
%% Frequency
    y=[];
for nidx=new_n_allruns
    y_ = zeros(30,size(cell2mat(allruns(nidx).curvature),2)); 
    for i = 1:size(cell2mat(allruns(nidx).curvature),2);
        y_(1:30,i) = allruns(nidx).curvature{i};
        %y = [y;allruns(nidx).curvature{i}];        
    end
    y = [y,y_];
    % 
%     for i = 1:size(cell2mat(allruns(nidx).curvature),2);
%         y = [y,allruns(nidx).curvature{i}];
%     end
    
%     y = [y;cell2mat(allruns(nidx).curvature)];
    %Fs =  params(1).time_res; %Hz
    Fs = params(1).time_res;
    
   
end

clear Y;
clear magnitudeY;

for i = 1:30

% y = y';
%     NFFT = length(y);
%     Y = fft(y(:),NFFT);

  
    NFFT = length(y);
    Y(i,:) = fft(y(i,:),NFFT);
    Y(i,:) = fftshift(Y(i,:));
    dF = Fs/NFFT;                      % hertz
    F = -Fs/2:dF:Fs/2-dF;           % hertz
    magnitudeY(i,:) = abs(Y(i,:));
    phaseY = unwrap(angle(Y(i,:)));
    %helperFrequencyAnalysisPlot1(F,magnitudeY,phaseY,NFFT)
    figure;
    plot(F, magnitudeY(i,:)); xlim([0 max(F)])
end

figure(90);
plot(F,sum(magnitudeY(1:30,:),1)); xlim([0 max(F)])

figure(91);
plot(F,sum(magnitudeY(1:3,:),1)); xlim([0 max(F)])

disp(['Number of tracks included:' num2str(length(new_n_allruns))])
